installed.packages()   #已经安装的包
.libPaths()           #路径查看
if(! require("ggplot2")) install.packages("ggplot2")
if(! require(" tidyverse ")) install.packages("tidyverse")
if(! require("ggplot2")) install.packages("ggplot2")
if(! require("ggplot2")) install.packages("ggplot2")
if(! require("ggplot2")) install.packages("ggplot2")
if(! require("ggplot2")) install.packages("ggplot2")

#基本命令包括表达式与赋值
#赋值  = 或<-  ，前者可传参或赋值，后者只可赋值
length(x=seq(1,10))   #传参给x，x为临时变量，环境中不可见
length(x<-seq(1,10))  #赋值给x,x为变量，可见

#大小写敏感，同一变量直接覆盖,变量名可以为中文，不可以下划线开头
中文=1
English=2
.abc=3
print(English)
print(中文,English,.abc)
print(ls())
cat(中文,English,.abc)   #马一诺
#print()只能打印一个变量
#cat()可以打印多个变量的组合【cat(“aaa”,var,”bbb”)】
#print(ls())打印工作区内所有的变量
#.开头为隐藏变量

#第2题
print("马一诺")
cat("马一诺")
print("monuo")
cat("monuo")



#条件语句
#if(A){B}
#if(A){B} else {c}
#ifelse(test,yes,no)

#循环语句
#for(i in a){语句}
i=1   #必须先定义,while语句中
while(i<10)
{print(i) 
  i=i+1}

#作业3-马一诺
num=sample(1:50,10,replace=TRUE)
cat("随机数组num[]为",num)    #输出num

#3.1
a=0   #小于等于10数的计数变量
i=1
for(i in num)   #循环法1
{if (i<=10)
{a=a+1}
}
b=10-a  #b为大于10的个数
print("             /n")
cat("大于10的有",b,"个")

#3.2
i=1    #重置循环计数变量
c=0    #小于40的计数变量
while(i < 11) #循环法2
{
if (num[i]<40)
{c=c+1}
i=i+1 
}
print("              /n")
cat("小于40的有",c,"个")



#3.3
#思路1:双层循环，
#外层循环9次，每次将剩下最大的数排在最前边；
#内层循环循环剩余数，找出其中最大数
num=sample(1:50,10,replace=TRUE)
cat("随机数组num[]为",num)    #输出num

for (m in (1:9))
{ 
  y=0
  for(n in (m:10))  #剩余的数
  {y=max(num[n],y)} #剩余的最大值
  x=m+which.max(num[m:10])-1 #数组最大值的索引
  num[x]=num[m]           #交换位置
  num[m]=y                
}
print(num)

#思路2 用sort函数直接排序做检验
#我本来想用“sort”作为数组名，后发现其是自带的排序函数，
#因此我改用“num”，并在学习sort函数后用其进行了结果验证
sort(num,decreasing = TRUE)

#思路三，调试完成后发现需要是条件语句，只能重新写max函数
for (m in (1:9))
{ 
  y=0
  for(n in (m:10))  #剩余的数
  {nummax=0
  if (nummax<num[m])
  {nummax=num[m]}
  } #条件语句求剩余的最大值,思路一中的max函数
  x=m+which.max(num[m:10])-1 #数组最大值的索引
  num[x]=num[m]           #交换位置
  num[m]=nummax                
}
print(num)

#4.22 实习三
#向量，一维数据；矩阵，二维数据；列表，数组  多维数据
#向量
A=c(1:10)     #小写c赋值
print(A)
B=seq(from=a,to=b,by=c) #从a到b,间隔为c
C=rep(0,5)         #重复5个0

#修改向量函数
A[-3]   #删除（用减号）第三个
A[-c(7,8)]  #同时删除第七，第8个
A[A<5]          #输出A向量里面值小于五的所有数
B=seq(from=1,to=10,by=2)
print(B[B<5])

#数值型向量函数常用的函数sum,max,sort，mean
#Inf -Inf  无穷数
#NA not available   缺失值--某个数值不可利用
#NaN not a number   缺失值--某个值不是一个数
#NULL  空值,不是缺失值

#字符型
#字符串处理函数：
#strsplit()分割,paste()及paste0()字符串连接
#length()向量长度；nchar()字符串长度

#逻辑型   TRUE、FALSE、NA

#一维数据--因子
#factor(x=character(),levels,labels=levels,exclude=NA,
#      ordered=is.ordered(x),nmax=NA)
#levels:用于指定因子可能的水平
#labels:用来指定水平的名字
#exclude:表示从向量x中剔除的水平值

#tapply()   基于因子的运算
gender=c('f','m','f')
age=c(10,11,12)
tapply(age,gender,mean)

#1.一维数据练习 马一诺
score=sample(85:100,10,replace=TRUE) #数值向量
group=c('富强','民主','和谐','民主','和谐','民主','民主','和谐','富强','富强')
group2=factor(group)
tapply(score,group,mean)
print(score)

factor(group,levels = '富强','民主','和谐',labels = 'A','B','C')
#将富强民主和谐换成A,B,C

#判断数据类型  is.double()
#转换数据类型  as.numeric()  转换成数值型

#二维--矩阵（里边数据类型一样）
mat=matrix(c(1:12),nrow,ncol=4)  #3行4列
#维数获取：dim()
#*  矩阵与标量×
#  %*%    矩阵之间×

#二维数据--数据框（里边可多种数据类型，但同一列数据类型相同）
dat=data.frame()
#$  数据框的每一列都有索引
#A$a  即数据框A的第a列
#rbind()  行合并，要求列相同    
#cbind()   列合并，要求行相同
#merge()  连接


#2.二维数据练习，马一诺，法1
score_R=sample(90:100,10,replace=TRUE) #数值向量2
score_python=sample(60:80,10,replace=TRUE)#数值向量3
d=data.frame(
  R分数=score_R,
  python分数=score_python,
  c语言分数=score,
  小组组别=group,
  stringsAsFactors=FALSE
)
knitr::kable(d)    #表格格式输出

#法2
data=cbind(score,score_R)
data2=cbind(data,score_python)
data3=cbind(data2,group)
knitr::kable(data3)    #表格格式输出

#法3
dataa=rbind(score,score_R)
datab=rbind(dataa,score_python)
datac=rbind(datab,group)
knitr::kable(datac)    #表格格式输出

#法2
data=cbind(score,score_R)
data2=cbind(data,score_python)
data3=cbind(data2,group2)
print(data3)    #表格格式输出

#函数
#形式： 变量名=function(变量列表) {函数体}

#输入输出

#路径
getwd()    #当前路径
setwd("")    #设置新路径，路径中的\要改成/
setwd
dir.create("路径/文件夹名")  #新建一个文件夹
list.files("路径")     #查询路径下的所有文件 
#3.2需要 可以专门读csv
#获取所有文件名、数据合并、删除多余信息


#read.table  read.csv  不能read.excel需要安装专用包   读
#write.table       写
#save.image(file="文件夹名")  把所有变量存进文件夹
#load()
#读文件前需要先确定文件类型和分隔符再进行读取。

#3文件的读入读出

#########
if(! require("openxlsx")) install.packages("openxlsx")
if(! require("readr")) install.packages("readr")  #安装readr包
library(readr)
setwd("C:/Users/34863/Desktop/data for class3")    #设置新路径
d2=read_csv("ADdata2.csv",locale=locale(encoding="GBK"))
d1=read_tsv("ADdata1.txt",col_names=FALSE)
d4=read_table2("ADdata4.txt")
if(! require("Rcpp")) install.packages("Rcpp")
library(openxlsx)
d3=read.xlsx("ADdata3.xlsx")
d6=read.csv("ADdata2.csv")
d5=read.delim("ADdata1.txt",header = TRUE,col_names = FALSE)

d8=data.frame(read.table("ADdata4.txt",header=TRUE))
names(d1)
names(d4)              #读取所有列索引
rownames(d8)           #读取所有行索引
names(d8)
names(d6)[1]=''
m1=merge(d5,d6,by.x='',by.y='')
m2=merge(d5,d6)
help(names)
###########   以上均为学习摸索的训练代码，下方为解题思路

# 3.1 读取并合并文件 马一诺
#① merge函数为主
#思路：依次读入四个文件后，找到共有列索引作为合并依据合并
#
library(openxlsx)            #掉包
setwd("C:/Users/34863/Desktop/data for class3")    #设置新路径
d1=read.delim("ADdata1.txt")
d2=read.csv("ADdata2.csv")
d3=read.xlsx("ADdata3.xlsx")
d4=read.table("ADdata4.txt",header=TRUE)
#读取列索引,可见列数和基因列索引各不相同，这是因为数据源首行首列为空导致
names(d1)
names(d2)
names(d3)
names(d4)
#思路，将基因名列添加数据框的一列，并将该列列索引统一为"X"
rowindex1=rownames(d1)      #读取行索引
rowindex4=rownames(d4)
d1new=cbind(d1,rowindex1)    #利用cbind函数将基因名添加至数据框
d4new=cbind(d4,rowindex4)    
names(d1new)[17]="X"         #统一基因名列索引
names(d4new)[17]="X"
names(d3)[1]="X"
#merge（）函数合并
data1=merge(d1new,d4new,by.x = "X",by.y = "X")
data2=merge(d3,d2,by.x = "X",by.y="X")
data3=merge(data1,data2,by.x="X",by.y="X")
#输出
write.csv(data3, file = "outcome1.csv")
write.table(data3, file = "outcome2.txt")
write.table(data3, file = "outcome3.xlsx")
#讨论：为了用R解决首行首列为空花费了不少时间，其实用excel等软件交互式的
#操作可以更便捷的添加字段，另外全部的整合过程excel也只需要导入数据，基因名排序，
#删除重复列三步即可。

#②cbind(),
#思路：直接按基因名列排序然后cbind，最后去除多余基因列
data4=cbind(d1new[order(d1new$X),],d4new[order(d4new$X),],d2[order(d2$X),]
            ,d3[order(d3$X),])    #这里要加，

#找出并去除重复列
names(data4)    #得到分别为17，34，35，52列
data5=data4[c(-34,-35,-52)]



###############################################################################################################################################
#实习课4，4/29日
#read.table()  header=TRUE,会自动把首行首列作为行名和列名
#3.1可以用（，rownames=null）完成
#read.table(check.names= TRUE) 会检查并自动修改列名
#存在唯一相同列名时，直接用merge(d1,d2) 就行
#不存在时可以采用列名不同的两列merge（d1,d2,by.X="AA",by.Y="BB"）即可
#甚至可以依照列名和某列来合并！
#merge（d1,d2,by.X="row.names",by.Y="BB"）即可

#数据查看：class 类型 dim 维度  head 快速查看前几行几列

#match 函数    
#match(x2,x1)  返回x2中所有元素在x1中的位置索引,不存在则返回NA
#x2  %in%  x1   返回x2每个元素是否在x1中



#gsub函数：
gsub("[a-zA-Z]+: ","",inf$tissue)   
#最后边是字符框与列，然后用中间的参数（null）替换首参
gsub("[ab]","",c("a","acashc","bbbba","abcba"))  #例子




#1数据整合
#1.1 马一诺 （上周是用order做的，比较便捷，这周采用循环依行添加的方法）
#数据用上次的方法读入并初步修改为d1new,d2,d3,d4new
row4=match(d1new$X,d4new$X)  
row2=match(d1new$X,d2$X)
row3=match(d1new$X,d3$X)
m1=d1new[1,]
m2=d2[row2[1],]
m3=d3[row3[1],]
m4=d4new[row4[1],]

qwer2=cbind(m1,m2,m3,m4)    #这里先用第一列试了一下，顺便获得首行添加结果

i=2
for(i in 2:100)
{
  qwer1=cbind(d1new[i,],d2[row2[i],],d3[row3[i],],d4new[row4[i],])
  qwer2=rbind(qwer2,qwer1)
}

#此时就完成了d1-d4的合并 
#接下来去需去除重复列即可
match("X",names(qwer2))     #????想用新学的match找到基因列
names(qwer2)          #想用新学的match，但只返回了一个值,只能用旧方法
qwer3=qwer2[c(-18,-35,-68)]




#list.files("路径"，pattern='AD")  只读含AD的  
 
#1.2需要 可以专门读csv  马一诺 在实习3的基础上改错与优化
#获取所有文件名、数据合并、删除多余信息
#为防止成字符型，最好不要rbind()，应转置后用merge


#字符操作：strsplit(flename,"_")，即将字符串按"_"分隔开

#获取所有文件名，并读取数据与合并
filename=list.files("C:/Users/34863/Desktop/data for class3/GSE67835")
setwd("C:/Users/34863/Desktop/data for class3/GSE67835")    #设置新路径
num_old=read.csv(filename[2],sep="\t",header = FALSE)
for(n in 3:139)                         #循环读入所有csv文件并连接
{num_new=read.csv(filename[n],sep="\t",header = FALSE)
num_old=cbind(num_old,num_new$V2)   #因为顺序是相同的
}
names(num_old)=(1:139)   #重命名列名

#读入txt文件
#先尝试直接读入看看前多少行不是分隔符
txt1=read.delim(filename[1])    #应跳过前34行（老师是跳35行，我觉得都行）
txt2=read.table(filename[1],skip=34,sep="\t",blank.lines.skip = TRUE,fill=TRUE)
#跳过空白行，填充null字段

sort(filename)==filename   #判断出基因数据是递增的
sort(txt2[2,])== txt2[2,]  #同样判断出txt的参数数据是递增的，
#因此可以直接合并，但老师说为了防止成字符型还是转置后用merge（）比较好

#提取目标行、转置、重命名
table2=txt2[c(2,10,11,12),-1]   #这里加-1确保为提取行，默认为提取列
table2=t(table2)  #转置
table2=data.frame(table2)
names(table2)=c("GSMnumber","tissue","celltype","age") #重命名

#利用gsub()函数去除冗余
table2$tissue=gsub("tissue: ","",table2$tissue)
table2$celltype=gsub("cell type: ","",table2$celltype) 
#注意表格中的cell type 是分开的
table2$age=gsub("age:","",table2$age)

#csv整合文件的重命名、转置
#之前用sort()检测过均为升序，一一对应
names(num_old)[2:139]=table2$GSMnumber  #简单很多，不再需要strsplit()
names(num_old)[1]="genename"
rownames(num_old)=num_old$genename
a=num_old[c(-1)]
table1=data.frame(t(num_old[c(-1)]))
table2
finaltable=merge(table2,table1,by.y="row.names",by.x="GSMnumber")

library(openxlsx)
write.xlsx(finaltable,"GSE.xlsx")

save(n3,num_old,txt2,table1,table2,finaltable, file="实践课4的GSE文件.RData")



#apply  将函数应用到数组维度
#apply(x, MARGIN, FUN, ...)
#其中，x为数据对象，MARGIN是维度的下标，FUN是由你指定的函 数，而...则包括了
#任何想传递给FUN的参数。在矩阵或数据框 中，MARGIN=1表示行，MARGIN=2表示列。

#2.马一诺
#思路，利用TRUE=1，FALSE=0的特性将问题转化为
#①布尔型判断每行每列是否为零
#②每行每列相加
x=function (h)   #定义计数函数x()
{sum(h==0)}

qwer3=qwer3[c(-17)]    #去除基因名列，防止字符型影响计算

rowmissing=apply(qwer3,1,x)      #每行缺失
columnmissing=apply(qwer3,2,x)   #每列缺失
columnmissing
rowmissing
#######################################################################################################################################

#######################################################################################################################################

#5.6   实习五
load

row4=match(d1new$X,d4new$X)  
row2=match(d1new$X,d2$X)
row3=match(d1new$X,d3$X)
m1=d1new[1,]
m2=d2[row2[1],]
m3=d3[row3[1],]
m4=d4new[row4[1],]


#read.table()  header=TRUE,会自动把首行首列作为行名和列名
#3.1可以用（，rownames=null）完成
#read.table(check.names= TRUE) 会检查并自动修改列名
#存在唯一相同列名时，直接用merge(d1,d2) 就行
#不存在时可以采用列名不同的两列merge（d1,d2,by.X="AA",by.Y="BB"）即可
#甚至可以依照列名和某列来合并！
#merge（d1,d2,by.X="row.names",by.Y="BB"）即可

#路径
getwd()    #当前路径
setwd("")    #设置新路径，路径中的\要改成/


data3=merge(data2,data1)
row1=match(data1$ID,data2$ID)
row1
row2=data1$ID-row1
m1=data1[1,]
m2=data2[row1[1],]
for(i in 1:100)
{
  qwer1=cbind(d1new[i,],d2[row2[i],],d3[row3[i],],d4new[row4[i],])
  qwer2=rbind(qwer2,qwer1)
}


data4=cbind(data1[order(data1$ID),],data2[order(data2$ID),])

data4=dplyr::bind_rows(data1[order(data1$ID),],data2[order(data2$ID),])
data4=dplyr::bind_rows(data1,data2)
library("plyr")

rownames(data2)=data2$ID
rownames(data1)=data1$ID
data5=plyr::rbind.fill(data2,data1)


install.packages("plyr")

data3=merge(data1,data2,all=TRUE)

#去重函数unique
#d5=unique(d5)


#数据处理
#cut(x,breaks=5)   将数值x分成五份，数据框的话要变成向量
#cut(unlist(x),breaks=5)   将行或列名去除，只保留数值


n3.1=unlist(n3[,1])
n3.1
#画直方图
hist(n3.1)
names(n3)[17]
ADdata=n3[c(-17)]
hist(n3.1,labels = TRUE)    #标出每列的数值


###########################################################
#1.1 马一诺
#ADdata
ADdata=n3[c(-17)]
ADdata2=unlist(ADdata)
ADdata3=ADdata2[!(ADdata2==0)]   #去除0值
#取对数
ADdata4=log2(ADdata3)
hist(ADdata4,label=TRUE,breaks=15,
     col=c("red","lightblue","green","orange"))

#GSE
GSE=unlist(table1)
GSE2=GSE[!(GSE==0)]   #去除0值
#取对数
GSE3=log2(GSE2)
hist(GSE3,label=TRUE,breaks=9,
     col=c("red","lightblue","green","orange"))

#######################################################

hist(ADdata3,label=TRUE,breaks=c(min(ADdata3),1e7,1e8,1e9,max(ADdata3)),
     col=c("red","lightblue","green","orange"))



include.lowest = FALSE  #去掉最小值
freq=F    #展示概率密度
col=c("red","lightblue")   #填充颜色


#1.2  马一诺
#选取P04004基因
#ADdata
P04004=log2(ADdata[7,])   #提取指定行,并取对数
P04004=unlist(P04004)
hist(P04004,label=TRUE,breaks=7,
     col=c("red","lightblue","green","orange"))

#选取A2M基因
#GSE
GSE=unlist(table1$A2M.)
A2M=log2(GSE[!(GSE==0)])   #去除0值并取对数
hist(A2M,label=TRUE,breaks=7,
     col=c("red","lightblue","green","orange"))






#PCA   主成分分析 （通过线性变换，实现降维）
#维度过多，纳入的无关因素过多时，往往会出现过拟合，假阳性的现象

#princomp()     observastion 需要大于 variable

prcomp

#标准化  

#ADdata先做转置，添加一列type看是asym,还是ad，分析asym之间的相似性
#PCA流程：去除无用列   rotation旋转和拉伸，维数不变，变为PC1-PCn
#选取PC1,PC2来画（这两维就基本反映了之前的四维）
#autoplot(col="type")  依据type 填成3色

#2.1  马一诺
if(! require("ggfortify")) install.packages("ggfortify")
#以下均为因4.0安装过需要重新装的包
if(! require("gtable")) install.packages("gtable")  
if(! require("colorspace")) install.packages("colorspace")
if(! require("labeling")) install.packages("labeling")
if(! require("digest")) install.packages("digest")
library("ggfortify")
ADdata2=t(ADdata)   #转置
#接下来提取asym、ad、ctl我有两种思路，一种是排序后分三批命名
#另一种是用讲过的gsub函数用null替换掉期中的数字。
#最后边是字符框与列，然后用中间的参数（null）替换首参
type=gsub("[0-9]","",rownames(ADdata2))     
#ADdata3=ADdata2[order(rownames(ADdata2)),]  这种方法也很简单
ADdata4=cbind(ADdata2,type)   #添加type列
ADdata.pca=prcomp(ADdata2,scale=T,rank=4,retx=T)
#标准化、返回score
autoplot(ADdata.pca,data = ADdata4, col= "type", size=2, loadings =T, 
         loadings.label = T, frame = T, 
         frame.type='norm', label = TRUE, label.size = 3) + 
  theme_classic()


#2.2 马一诺
GSE.pca=prcomp(finaltable[,5:22092],scale=F,rank=4,retx=T)
gc() #下面的处理将占用很大的内存，在运行之前先删去不必要的
autoplot(GSE.pca,data = finaltable, col= "celltype", size=2, loadings =F, 
         loadings.label = T, frame = T, 
         frame.type='norm', label = TRUE, label.size = 3) + 
  theme_classic()
#因为太多了，我将loadings改为了FALSE，因为具有全0列，我将scale改为了FALSE



memory.size(T)
memory.size(F)
memory.limit()
#######################################################################################################################################
#5.13   实习六

GSE=finaltable
load("class6_volcano.RData")
save(ADdata,file = "ADdata.RData")
save(GSE,file="GSE.RData")
load("GSE.RData")



mycol=c("#FF0000",#FFFF00","#008B8B",”#7FFFD4″,”#FFFAFA”,”#0000FF”,
          "#8A2BE2","#A52A2A","#000000","#7FFF00","#80000040","#FF7F50",
          "#6495ED","#DC143C","#00FFFF","#B8860B","#A9A9A9","#006400","#FFDAB9",
          "#8B008B","#FF00FF","#483D8B","#2F4F4F","#D2B48C"
)
save(mycol,file="mycol.RData")
#6.1  马一诺
#层次聚类
#dist函数
#https://blog.csdn.net/woodcorpse/article/details/84899008

dist_ADdata=dist(t(ADdata),method="manhattan")  
#曼哈顿法求距离，每一行之间的距离
dist_ADdata=hclust(dist_ADdata,method="average")
#hculst,根据距离分类,有average,centroid,median多种
plot(dist_ADdata,hang=-1,cex=.8,col=c("red","lightblue","green","orange"))  
#hang=-1会对齐最下方  ，cex.8是点的大小




#火山图绘制    马一诺
load("class6_volcano.RData")
#https://www.jianshu.com/p/e2828ef9c7e5
#加载包
library(ggplot2)
#abs 绝对值函数
foldchange=2^prostat$FC    #log2的反运算
#分三种情况填入新列
prostat$threshold = (ifelse(prostat$P < 0.05 & abs(foldchange) > 1.2, 'UP',
                                  ifelse(prostat$P < 0.05 &foldchange<= 1/1.2 ,
                                         'Down','NoSignifi')))
#采用ifelse的嵌套添加新列，标识基因是否特异表达(刚学的，觉得挺牛)

prostat$ID2=(ifelse(prostat$threshold != "NoSignifi",prostat$ID,''))
#添加ID2列，单独记录上行表达和下行表达的基因名，

jpeg("8305190706_monuo_火山图.jpg")    #保存
table(prostat$threshold)  
#table查看发现是down、nosignifi、up的顺序，指示后来颜色的命名
ggplot(prostat,aes(x=FC,y=-log10(P),color=threshold))+
  geom_point()+
  scale_color_manual(values=c("green","gray","orange"))+#确定点的颜色
  xlab('log2 (FoldChange)')+#修改x轴名称
  geom_vline(xintercept=c(-log2(1.2),log2(1.2)),lty=3,col="black",lwd=1) +
  #添加竖线|FoldChange|>1.2
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=1)+
   #添加横线P<0.05 
  geom_text_repel(aes(label=prostat$ID2),color="black",size=3)
#添加基因标签
dev.off()


prostat$threshold
1!=1


install.packages("ggrepel")
library(ggrepel)
  geom_text_repel(aes(label=FC),color="black",size=3)
#添加基因标签
a


load("mycol.RData")
1/1.2
pdf("name.pdf")
#画图语句
dev.off

#alpha 透明度
#作图
ggplot(prostat,aes(x=FC,y=-log10(P)))+
  geom_point()
ggplot(Dat,aes(x=log2FoldChange,y=-log10(padj),color=threshold))+
  geom_point()+
  scale_color_manual(values=c("#DC143C","#00008B","#808080"))+#确定点的颜色
  geom_text_repel(
    data = Dat[Dat$padj<0.05&abs(Dat$log2FoldChange)>1,],
    aes(label = Gene),
    size = 3,
    segment.color = "black", show.legend = FALSE )+#添加关注的点的基因名
  theme_bw()+#修改图片背景
  theme(
    legend.title = element_blank()#不显示图例标题
  )+
  ylab('-log10 (p-adj)')+#修改y轴名称
  xlab('log2 (FoldChange)')+#修改x轴名称
  geom_vline(xintercept=c(-1,1),lty=3,col="black",lwd=0.5) +#添加横线|FoldChange|>2
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=0.5)#添加竖线padj<0.05




#mean+-SD
ad=ifelse(type="ad",ADdata[2,],null )
colnames(ADdata)

type=gsub("[0-9]","",names(ADdata))
table(type)
table(names(ADdata))
b=a[order(rownames(a)),]
a=t(ADdata[2,])
b=b[-40:-54]
#######################################################################################################################################

#######################################################################################################################################
#5.20   实习7  马一诺
load("ADdata.RData")
load("class6_volcano.RData")

#显著性检验：T检验、方差检验、ANOVA、线性回归、
#wilcoxon秩和检验  （处理非正态分布数据）
#limma 基于线性回归模型的Empirical Bayes参数估计
#P值：是一样的概率。所以越小差异越显著

#t检验
t.test(x,y=NULL,alternative=c("two.sided","less","greater"),
       mu=0,paired=FALSE,var.equal=FALSE,conf.level=0.95,...)
#只有X单样本，有X、Y二样本t检验
#alternative=c("two.sided","less","greater")用于指定所求置信区间的类型
#alternative="two.sided" 是缺省值，表示求置信区间

#pairwise.t.test(pool.sd=FALSE,var.equal=TRUE) 多变量t检验
#正态分布检验    
#shapiro.test()
#方差齐性检验
#x=aov(GLU~sex,diabets)


#实习七 1.  马一诺


#要3类都至少有3个非0的
#1 对数化  1.1去0，把0值变成NA  1.2对数化

#2  提取1个蛋白的数据，提取分组信息
#3  判定是否满足单组非0值>3
#4   方差分析或输出NA
#5   如何迭代到循环里

#1.  马一诺
ADdata2=t(ADdata)   #转置
ADdata2[ADdata2==0]=NA  #去0,将0转换成NA
ADdata2=log2(ADdata2)   #取对数
type=gsub("[0-9]","",rownames(ADdata2))   
c=names(as.data.frame(ADdata2))  #先定义一个100的，用来存结果
#接下来我的思路和老师略有不同
#只要去除type列中对应NA的部分再用table函数计数即可，
#不用分三组分别提取对应值在计数
i=1
for(i in (1:100))
  {
ADdata3=cbind(ADdata2[,i],type)     #提取type和第i列
ADdata3[,2]=ifelse(is.na(ADdata3[,1])==T,NA,ADdata3[,2])
#ifelse语句对应NA处
b=table(ADdata3[,2])>3  #存储三种type的类型
{if  (all(b)==T)     #这里学习了all和any
{result=oneway.test(ADdata2[,i]~type)
c[i]=result$p.value
}
  else
  {
    c[i]=NA
  } } 
}




c



#if(A){B} else {c}

result$p.value
class(result)













#老师讲的地方
a=as.vector(ADdata2[,1])   #变成向量

posad=(type=="ad")  #ad是哪些
ad=a[posad]            #返回ad的数
#同理得到三组
ad2=!is.na(ad)   #哪些不是NA
adsum=sum(ad2)   #计数不是NA的数量

#if  (adsum>3  &  )  #(&是且，|是或)    
{result=oneway.test(a~type)}  #算p值

d=as.data.frame(result)



ADdata2=t(ADdata)   #转置
rownames(ADdata2)=gsub("[0-9]","",rownames(ADdata2))    

ADdata3=cbind(ADdata2,type)   #添加type列

data1=unlist(ADdata2[,1])
data1=log2(data1[!(data1==0)])   #去除0值并取对数
for(i in 100)   #循环法1
{data=unlist(ADdata2[,i])
data=log2(data[!(data==0)]) 
{}
if (length(data[i])<4)
{
  p=NA
}
}


class(data1)
data1[1]
rownames(ADdata2)

length(data1)


#利用矩阵乘法快速算出有多少0
data1=ADdata
detectnum=(t(data1)!=0)%*%rep(1,nrow(data1))
#pos看0都在哪
pos=which(detectnum==0);
data1=data1[,-pos]
data1$P01008

data=GSE[,-1:-4]
detectnum=(t(data)!=0)%*%rep(1,nrow(data))
#pos看0都在哪
pos=which(detectnum==0);
data=data[,-pos]


#相关性分析
#富集分析


#ebrichGO()  富集
#ont="" ,类型 ，universe 最好设成背景基因 
#barplot()   柱状图展示
#dotplot(ena2)  气泡图
#emapplot()       基因关联图
#https://blog.csdn.net/weixin_43569478/article/details/83744242

if(! require("clusterProfiler")) install.packages("clusterProflier")
if(! require(" tidyverse ")) install.packages("tidyverse")

source("https://bioconductor.org/biocLite.R")
biocLite("DOSE") 
biocLite("topGO")
biocLite("clusterProfiler")   
biocLite("pathview")
library(DOSE)
library(org.Hs.eg.db)       #专指人类的go数据库
library(topGO)
library(clusterProfiler)
library(pathview)




if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(version = "3.14")
BiocManager::install("clusterProfiler")
n
library("clusterProfiler")


options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos" = c(CRAN="http://mirrors.cloud.tencent.com/CRAN/")) 
options(download.file.method = 'libcurl')
options(url.method='libcurl')
BiocManager::install("GO.db",ask = F,update = F)
BiocManager::install("org.Hs.eg.db",ask = F,update = F)
#实习7  第2、3题  马一诺
library(org.Hs.eg.db)
library(clusterProfiler)
library(ggplot2)
#加载数据
load("class6_volcano.RData")
#数据预处理，筛选p值小于0.05的蛋白
P = prostat$P
i= P < 0.05
pro_data = prostat[i,]

allID = pro_data$ID[!is.na(pro_data$ID)]
enall = enrichGO(allID,"org.Hs.eg.db",keyType = "SYMBOL",ont = "ALL",pvalueCutoff = 0.05,
                  pAdjustMethod = "BH",minGSSize = 10,maxGSSize = 500)

#条形图展示结果
barplot(enall, showCategory = 20,drop = T)
#气泡图展示结果
dotplot(enall, showCategory = 30)







update.packages()



if(! require("Biobase")) install.packages("Biobase")

if(! require("Biobase")) install.packages("Biobase")

if(! require("Biobase")) install.packages("Biobase")

if(! require("Biobase")) install.packages("Biobase")

install.packages("Biobase")
BiocManager::install("Biobase")







if(! require("psych")) install.packages("psych")

load("pro_data.RData")
load("add_file.RData")
load("BLSA_data.RData")

library(psych)
add_file.des <- describe(add_file,IQR=TRUE)
#blsa_pro:
#直方图
hist(unlist(blsa_pro[,-1:-3]),rainbow(15))
#PCA
library(ggfortify)
blsa_pro.t <- as.data.frame(t(blsa_pro[,-1:-3]))
label <- rownames(blsa_pro.t)
blsa_pro.t[,"label"] <- label
blsa_pro.t$label <- gsub("_[0-9]+","",blsa_pro.t$label)
blsa.pca<- prcomp(blsa_pro.t[,-ncol(blsa_pro.t)],scale=T)
autoplot(blsa.pca,data=blsa_pro.t,col = 'label',size=0.5,frame=T,
         frame.type='norm')+theme_classic()





#火山图
load("ad_ct.RData")
table(pro_data$group)
add_file=add_file[-35,]
data1=cbind(add_file,pro_data)


load("class6_volcano.RData")
#https://www.jianshu.com/p/e2828ef9c7e5
#加载包
library(ggplot2)
library(ggrepel)
#abs 绝对值函数
foldchange=2^ad_ct$log2FC    #log2的反运算
#分三种情况填入新列
ad_ct$threshold = (ifelse(ad_ct$P < 0.05 & foldchange > 1.5, 'UP',
                            ifelse(ad_ct$P < 0.05 &foldchange<= 1/1.5 ,
                                   'Down','NoSignifi')))
#采用ifelse的嵌套添加新列，标识基因是否特异表达(刚学的，觉得挺牛)

ad_ct$ID2=(ifelse(ad_ct$threshold != "NoSignifi",ad_ct$ID,''))
#添加ID2列，单独记录上行表达和下行表达的基因名，


table(ad_ct$threshold)  
#table查看发现是down、nosignifi、up的顺序，指示后来颜色的命名
ggplot(ad_ct,aes(x=log2FC,y=-log10(P),color=threshold))+
  geom_point()+
  scale_color_manual(values=c("green","gray","orange"))+#确定点的颜色
  xlab('log2 (FoldChange)')+#修改x轴名称
  geom_vline(xintercept=c(-log2(1.5),log2(1.5)),lty=3,col="black",lwd=1) +
  #添加竖线|FoldChange|>1.5
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=1)+
  #添加横线P<0.05 
 geom_text_repel(aes(label=ad_ct$ID2),color="black",size=3)
#添加基因标签



#ADdata2[==0]=NA  #去0,将0转换成NA


label=ad_ct$threshold
up=label!="NoSignifi"
ad_change_protein <- ad_ct[up,]
save()

save(ad_change_protein,file = "ad_change_protein.RData")

install.packages("tidyverse")
library(tidyverse)

tad <- t(pro_data)
rown <- rownames(tad)
tad$ID <- rownames(tad)

#  树状图
dist_ADdata=dist(pro_data,method="manhattan")  
#曼哈顿法求距离，每一行之间的距离
dist_ADdata=hclust(dist_ADdata,method="average")
#hculst,根据距离分类,有average,centroid,median多种
plot(dist_ADdata,hang=-1,cex=.8,col=c("red","lightblue","green","orange"))  
#hang=-1会对齐最下方  ，cex.8是点的大小

#######################################################################################################################################